function [ Prop,  RT, RT_STD ] = GaussSimulateRT( d, c, mu, sigma)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

Prop = zeros(length(c),2);
RT = zeros(length(c),2);
RT_STD = zeros(length(c),2);

for i = 1 : length(c)
    [Prop(i,:),  RT(i,:), RT_STD(i,:)] = GaussRandomDots(d, c(i), mu, sigma);
end

end

function [Prop, RT, RT_STD] =  GaussRandomDots(d, c, mu, sigma)
numTrials = 10000;
P = 0;
N = 0;
iter = 0;
R_P = 0;
R_N = 0;
STD_P = 0;
STD_N = 0;
nT = size(d,1);
thrR = d(:,2)';
thrL = d(:,1)';
while iter < numTrials
    iter = iter + 1;
    o = mu * c + randn(1,nT) * sigma;
    %o(1) = 0;
    o = cumsum(o) ./ (1:nT);
    tR = find(o >= thrR);
    tL = find(o <= thrL);
    if isempty(tR)
        tR(1) = nT;
    end

    if isempty(tL)
        tL(1) = nT;
    end
    
    if(tR(1) > 1)
        rB = tR(1) - 1 + (o(tR(1)-1) - thrR(tR(1)-1)) / (thrR(tR(1)) - thrR(tR(1)-1) - o(tR(1)) + o(tR(1)-1) );
    else
        rB = 1;
    end
    if tL(1) > 1
        lB = tL(1) - 1 + (o(tL(1)-1) - thrL(tL(1)-1)) / (thrL(tL(1)) - thrL(tL(1)-1) - o(tL(1)) + o(tL(1)-1));
    else
        lB = 1;
    end
    
    if tR(1) < tL(1)
        P = P + 1;
        R_P = R_P + rB;
        STD_P = STD_P + rB * rB;
    end    
    
    if tL(1) < tR(1)
        N = N + 1;
        R_N = R_N + lB;
        STD_N = STD_N + lB * lB;
    end
end
R_P = R_P / P;
R_N = R_N / N;
STD_P = sqrt((STD_P - R_P * R_P * P) / (P - 1));
STD_N = sqrt((STD_N - R_N * R_N * N) / (N - 1));
Prop = [P N];
RT = [R_P R_N];
RT_STD = [STD_P STD_N];

end